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We perform cosmological N-body simulations of the Dvali-Gabadadze-Porrati braneworld model, 
by solving the full non-linear equations of motion for the scalar degree of freedom in this model, 
the brane bending mode. While coupling universally to matter, the brane-bending mode has self- 
interactions that become important as soon as the density field becomes non-linear. These self- 
', interactions lead to a suppression of the field in high-density environments, and restore gravity to 

■ General Relativity. The code uses a multi-grid relaxation scheme to solve the non-linear field equa- 

' tion in the quasi-static approximation. We perform simulations of a flat self-accelerating DGP model 

, without cosmological constant. However, the type of non-linear interactions of the brane-bending 

mode, which are the focus of this study, are generic to a wide class of braneworld cosmologies. The 
results of the DGP simulations are compared with standard gravity simulations assuming the same 
expansion history, and with DGP simulations using the linearized equation for the brane bending 
mode. This allows us to isolate the effects of the non-linear self-couplings of the field which are 
• noticeable already on quasi-linear scales. We present results on the matter power spectrum and 

the halo mass function, and discuss the behavior of the brane bending mode within cosmological 
structure formation. We find that, independently of CMB constraints, the self-accelerating DGP 
model is strongly constrained by current weak lensing and cluster abundance measurements. 
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CLc I. INTRODUCTION 

O 

' The observed acceleration of the universe [1-5] poses a fascinating challenge to physicists and cosmologists: it 
C/3 , points towards new physics at an unusually low energy scale (~ 10~^eV), or large length scale (~ lO'^Mpc). For 
this reason, no natural explanation appears to exist up to now. Many attempts have been made to go beyond the 
minimal explanation, a cosmological constant or vacuum energy. They can be broadly classified into two categories: 
CnI acceleration is due to an additional, smooth stress-energy component with negative pressure (dark energy [6]); or 
I it is caused by gravity itself, through modifications to General Relativity (GR) on cosmological scales. While there 
are strong constraints on deviations from GR in the Solar System [7], gravity is remarkably weakly constrained on 
cosmological length scales. This provides an independent motivation for studying modified gravity models in the 
context of cosmology. 

Smooth dark energy models have enough freedom to reproduce essentially any background expansion history of 
the universe. In order to distinguish modified gravity from the smooth dark energy scenario, it is thus necessary 
to study the growth of cosmological structure. A wealth of observables can be used for this purpose [8, 9, 9-18], 
[ e.g. weak lensing, galaxy-CMB cross-correlation, galaxy cluster abundances, and many more. However, almost all of 
these measurements are affected by non-linearities in the matter density field, or have most of their information on 
non-linear scales. All viable modified gravity models include a non-linear mechanism to restore General Relativity in 
high-density environments, which is necessary in order to satisfy Solar System constraints. This mechanism has to be 
^ • taken into account in order to make accurate predictions for observables on non-linear scales in modified gravity. For 
5^ , one representative of f{R) modified gravity, [19, 20] showed that the effects of this non-linear chameleon mechanism 
are significant for models that satisfy Solar System constraints. 

One of the most popular modified gravity scenarios is the Dvali-Gabadadze-Porrati (DGP) model [21]. In this 
model, matter and radiation are confined to a 4-dimensional brane in a 5-dimensional bulk spacetime. While gravity 
propagates in five dimensions on the largest scales, it becomes 4-dimensional below a certain crossover scale re- On 
scales smaller than rc, and when gravity is weak, DGP gravity can be described as an effective four-dimensional 
scalar-tensor theory [22]. The scalar field, referred to as brane bending mode, is massless but has non-linear derivative 
interactions which suppress the field in high-density environments, restoring GR locally and allowing DGP gravity to 
pass Solar System constraints. 

Applied to homogeneous and isotropic cosmology, the DGP model allows for two solutions [23]. In one branch of the 
theory, the self- accelerating branch, the effective weakening of gravity on large scales leads to an accelerated late-time 
expansion of the brane, without any cosmological constant. The scalar brane-bending mode mediates a repulsive force, 
leading to a smaller effective gravitational constant on large scales. In the other solution, the normal branch, there 
is no accleration, and the brane-bending mode mediates an attractive force. In linear perturbation theory around its 
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dc Sitter limit, the brane bending mode in the self- accelerating DGP model has been shown to have the wrong sign 
in the kinetic term ("ghost") [22, 24], suggesting an instability, although the situation in the non- linear case is not 
clear [25, 26]. The normal branch is free of the ghost. 

Linear perturbation theory around a cosmological background in DGP has been derived in, e.g. [17, 27, 28]. In 
this approach, it is possible to predict the expansion history, CMB anisotropics, and the ISW effect in DGP [29]. 
Using this information, [30] have shown that the self-accelerating DGP model without cosmological constant, which 
has only one free parameter, rc, is disfavored at the ~ 4cr level by current CMB and Supernova data. This is due to 
the earlier onset of acceleration in DGP, and the additional suppression of growth through the brane-bending mode. 

While the self-accclcrating branch of DGP is thus challenged by theory and observations, much work is being 
done on extending the DGP model to higher co-dimensions, in the context of degravitation [31, 32]. These higher- 
dimensional models are expected to bring the expansion history close to that of ACDM, while exhibiting a similar 
effective scalar-tensor regime as in the original DGP model (see also [33]). The form of the non-linear interactions of 
the brane-bending mode are also generic to these generalized models (see also [34]). Hence, it will be straightforward 
to perform simulations of these models once the expansion history and evolution of linear perturbations is worked 
out. In addition, the standard DGP model is able to satisfy the cosmological constraints if a cosmological constant 
(brane tension) is included [35]. 

So far, studies of the DGP model in the context of cosmology have mostly dealt with the linearized theory, valid on 
large scales. However, as soon as perturbations in the matter density become of order unity, the non-linear interactions 
of the brane bending mode become important (e.g., [36]). Hence, in order to make predictions for observables in the 
non-linear regime, the full non-linear equations of motion must be solved in conjunction with the evolution of the 
matter perturbations. Recently, perturbative approaches have been developed to extend the predictions into the 
quasi-linear regime [37]. Also, Khoury and Wyman [38] have used a spherically-symmetric approximation of the 
brane-bending mode interactions in their N-body simulation of braneworld models. In this paper, we present the 
results of an N-body simulation of the self-accelerating DGP model, which self-consistently solves the full non-linear 
equation for the brane-bending mode and its effect on the motion of particles (we compare our results with the 
approximation used in [38] in Section VIC). We solve the equation of motion in the quasi-static approximation, 
neglecting time derivatives with respect to spatial derivatives, as is usually done in N-body simulations (Section IV C 
presents a consistency test of this approximation). 

We present results on the matter power spectrum, the halo mass function (abundance of dark matter halos), and the 
behavior of the brane bending mode in the cosmological context. These results can be used to strengthen constraints 
on the self- accelerating DGP model significantly, e.g., through observations of weak Icnsing and the abimdance of 
galaxy clusters. More generally, the simulations presented here can serve as starting point for building a model of 
non-linear structure formation in braneworld scenarios, and for benchmarking perturbative approaches [37] . 

In Section II, we describe the DGP model, and present the relevant equations and analytical test cases. Section III 
presents the code, which is benchmarked in Section IV. The cosmological simulations are described in Section V. 
Results and their impact on cosmological constraints on the DGP model are presented in Section VI and Section VII. 
We conclude in Section VIII. 



II. DGP MODEL 
A. Background evolution 

The Dvali-Gabadadze-Porrati model [21] consists of a spatially three-dimensional brane in a 4-1-1 dimensional 
Minkowski bulk. Matter and all interactions except gravity are confined to the brane. The gravitational action 
consists of the five- dimensional Einstein-Hilbert action plus a term which leads to the 4D gravity limit on small 
scales: 

^grav = -J^:^ J d'xV^R'-'^ + ^boundary - J^;^ J d^Xy/^ (R^"^ + . (2.1) 

Here, X, g^^^ stand for the bulk coordinates and metric, while x, g^^^ are the induced coordinates and metric on the 
brane, and R'^^\ i?*-^' denote the corresponding Ricci scalars. In the following, we will drop the notation where no 
confusion can arise. The boundary term is added to the action in order to ensure that variation with respect to g^^^ 
leads to the correct five-dimensional Einstein equations (e.g., [39]). 

The two gravitational constants, or Planck masses Mpi(5), Mpi(4) appearing in Eq. (2.1) can be related via a length 
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scale, the crossover scale: 

IG^_IMI,,^ (2.2) 



2GW 2M3^(,/ 

On scales above Tc, gravity becomes five-dimensional, with forces falling off as l/r^. Below r^ gravity is four- 
dimensional, but not Einsteinian gravity, a point to which we return below. 

Since all matter is thought as confined to the brane, the five- dimensional metric has to obey non-trivial junction 
conditions over the brane [23]. Assuming an empty Minkowski bulk, a spatially flat brane. and a homogeneous and 
isotropic matter distribution on the brane, the junction conditions lead to the following analogue of the Friedmann 
equation for the scale factor of the induced metric on the brane: 

Tc 3 a 

The junction conditions leave two possible branches of the theory, determined by the sign on the left hand side of 
Eq. (2.3). In the following, we focus on the branch with the '— ' sign, which asymptotes to a late-time de Sitter- 
Universe, H = l/fc = const., and is correspondingly called the accelerating branch. Then, in the matter-dominated 
epoch, neglecting radiation, and again assuming no cosmological constant or curvature, Eq. (2.3) can be rewritten as: 

H{a) 



E{a) = — ^ = + V^mCi-^ + ^rc, where (2.4) 

-no 

o = 1 o =^ 

and pmQ is the average matter density today. This expansion history is clearly different from ACDM, and corresponds 
to an effective dark energy with Wgs ~l/2 in the matter-dominated era at high redshifts. 



B. Cosmological perturbations and brane bending mode 

The propagation of light and particles on the DGP brane is completely determined by the perturbed 4D Friedmann- 
Robertson- Walker metric: 

ds"^ = -[1 + 2^{x,t)]dt^ + [1 + 2$(x,t)]a2(i)dx2 (2.5) 

However, in order to determine the evolution of the metric potentials on the brane, it is necessary to solve the full 5D 
Einstein equations [17, 36]. An additional scalar degree of freedom associated with local displacements of the brane 
appears, the so-called brane- bending mode ip which couples to matter. In our convention, ip is dimensionless, instead 
of being scaled to the Planck mass Mpi (4) . 

In the decoupling limit of DGP [22], when setting gravitational interactions to 0, the self- interactions of the (f field 
remain constant. Hence, while perturbations of the metric higher than linear order can be neglected for cosmological 
studies, it is crucial to consider the self-interactions in the brane-bending mode. 

In the quasi-static regime, for scales k ^ _ffo, i time derivatives can be neglected with respect to spatial 
derivatives, and the equation for the brane-bending mode reads (e.g., [36]): 

vV + ^[(vV)' - (V,V,-^)(V^VV)] = ^^^P. (2-6) 

where the derivatives are with respect to comoving coordinates x, Sp = pm — Pm is the matter density perturbation, 
and the function (3{a) (for the accelerating branch) is given by: 

/3(a) = l-2i^(a)..(l + Jg^). (2.7) 

Note that /3 is always negative in the self-accelerating branch. In Section IV C, we show results of a consistency test 
of the quasi-static assumption used in Eq. (2.6). 

The brane-bending mode influences the dynamics of particles through the dynamical potential which, assuming 
the same boundary conditions for ^ and ip, is given by: 

* = *iv + ^V', (2.8) 
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where the Newtonian potential satisfies the usual Poisson equation: 

V^^iv =47rGa^^p. (2.9) 

In contrast, the propagation of photons, determined by the lensing potential = ($ — ^')/2, is not directly affected 
by if. 

C. Weak brane regime 

If the gradient of the field is small, i.e. for small gravitational accelerations, Eq. (2.6) can be linearized, yielding 
a standard Poisson equation: 

„2 SwGa^, 2 ^ , 

Wc will refer to a cosmology with this linearized equation as the linearized DGP m,odel. In this regime, also called 
weak- brane phase, ip = (pL becomes proportional to the Newtonian potential ^Fat, corresponding to a constant rescaling 
of the gravitational constant through Eq. (2.8): 

Gn ^ Geff =(l + ^ Gn. (2.11) 

Substituting the linear solution into Eq. (2.6), we obtain a rough estimate for the overdensity 5 = 5p/pm at which 
the non-linear interactions become important: 

For a self-accelerating DGP model that leads to cosmic acceleration today, Tc ^ , (3{a = 1) ~ —1, and the prefactor 
in Eq. (2.12) is of order unity today. Hence, the self-coupling of the brane bending mode </? becomes important as 
soon as the matter density field becomes non- linear, (5^1. 

D. Vainshtein effect 

In general, Eq. (2.6) is difficult to solve in full generality due to the non-linearity in the derivative terms. However, 
an instructive test case, a spherically symmetric matter distribution, can be solved analytically. In the spherically 
symmetric case, Eq. (2.6) becomes [36]: 



f <P 2 d \ , 1 StxGo? , 



E ^ 2rl 



dr' / dif^ 







For simplicity, we assume a spherical mass M of radius R with uniform density, and set a = 1. Then, we can integrate 
Eq. (2.13) once and obtain the gravitational acceleration in DGP: 

9 = 9N + l^=9N[l + A{r)], (2.14) 



where ^jv is the Newtonian acceleration of the spherical mass, and: 

- 1 j , r>R 

(2.15) 
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Here, r* denotes a characteristic scale of the solution, the Vainshtein radius: 

and = 2GM is the gravitational radius of the mass. For very large distances, r ^ r*, A(r) approaches the constant 
1/(3/3), which exactly matches the linear solution, Eq. (2.11). Note also that by substituting 6p ~ M/r^ in Eq. (2.12), 
we see that the non-linearity criterion is directly proportional to (r*/r)^. 

In the opposite limit, r <C r*, 99 becomes suppressed with respect to the linear solution, and A(r) approaches 
the small constant 2/(3/3)(i?/r,)''/^ inside the mass (assuming R ^ r,). This Vainshtein effect [40, 41] amounts to 
restoring GR in deep potential wells, while far away from the mass, gravity is in the scalar-tensor regime. 

Fig. 1 shows the analytical solution for the ip field and the deviation from the Newtonian acceleration, as well as 
the code results which agree well with the analytical solution in the regime of validity (see Section IV A). 



E. Plane wave solution 



Another exact solution to Eq. (2.6) exists: for a simple plane wave density perturbation, Sp{x) = A exp(ik • x). In 
this case, the two non-linear terms exactly cancel, and we are left with the linear solution, Eq. (2.10). In Appendix B 

this is used to test our code for perturbations of different wavelengths. Also, the plane wave and spherically symmetric 
solutions can be considered as two limiting cases for understanding the (p field behavior in the cosmological context. 
We will return to this point in Section VI C. 



III. CODE IMPLEMENTATION 



We use an extension of the code first described in [42] for our simulations of DGP gravity. The code is a standard 
particle-mesh N-body code [43, 44] assuming collisionless dark matter only, with a fixed grid of size and periodic 
boundary conditions, and uses second-order accurate leapfrog integration for the particle propagation. The major 
addition is a relaxation solver for the non-linear equation of the brane bending mode (f, Eq. (2.6). The simulation 
proceeds by performing steps Aa in the scale factor a. 

Except for the different expansion history, the equations for particle propagation are identical to those in standard 
CDM simulations, and are expressed in terms of comoving coordinates x, p: 

'^"'^ ' (3.1) 



da a^ E{a) 
dpi 1 



da aE{a) 



V*(xi). (3.2) 



Here, and pj are the position and momentum of particle i in code units (see Appendix A), and aE{a) is taken 
from Eq. (2.4). The main modification enters in the determination of ^(x). One time step of the N-body simulation 
proceeds as follows: 



• The density field on a fixed grid is assigned using the particle positions at scale factor a via the cloud-in-cell 

Om ''cell- 



method, i.e. each particle corresponds to a uniform density cube of side length rceii and mass M = p^r^ 



• The Newtonian potential ^'tv is solved for using a Fast Fourier Transform; Eq. (2.6) is solved via relaxation, and 
the overall dynamical potential \l/ is assembled according to Eq. (2.8). We also run linearized DGP simulations, 
where the potential is simply given by ^' = [1 -|- l/(3/3(a))] ^'jv- 



Using the same cloud-in-cell interpolation as for particles, the gradient of 4" is calculated at each particle position 
and is used to update the particle momenta from a — Aa/2 to o -|- Aa/2. 

Using the momenta at a -|- Aa/2, the particle positions are updated from a to a -|- Aa, and the process starts 

from the beginning. 



Thus, in each time step we have to solve for the brane-bending mode (p which contributes to the dynamical potential 
via Eq. (2.8). In other words, given the density field, we have to solve the non-linear elliptical differential equation 
Eq. (2.6). We use a Gauss-Seidel relaxation scheme together with the Newton method. A crucial tool is the multigrid 
[45-47] , a hierarchy of coarser grids which are essential to speeding up the convergence [42] : as the relaxation scheme. 
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FIG. 1: Left panel: (p{r) measured in simulations with a 512^ grid for a spherical mass for different values of /3, where r is the 
distance from the center of the mass. R = 10 grid cells, and rc — 800 R is held fixed (see text). The solid lines show the full 
analytical solution for the spherical mass, while dotted lines show the linearized solution (Section II C). Right panel: Relative 
deviation of the acceleration measured in the spherical mass solution from the analytical Newtonian value vs. r, for the same 
parameters as in the left panel. Again, solid lines show the full exact solution, while dotted lines show the linearized solution. 



which operates locally, gets rid of small-scale errors very efficiently, the coarser grids quickly reduce the long-wavelength 
error modes which are hard to eliminate on the fine grid alone. For details of the implementation, see Appendix A. 

While the ip field is in general well-behaved, the convergence properties become worse for strongly inhomogeneous 
density fields, because the non-linearities in Eq. (2.6) arc in the derivatives of the ip field. In order to reach the desired 
convergence in our cosmological simulations, it is necessary to smooth the density field entering the r.h.s. of Eq. (2.6) 
with a Gaussian kernel cx exp(— r^/2r^), with set to the size of a grid cell. This smoothes over the noise in the 
density field due to the discreteness of particles. 



IV. CODE TESTS 



In this section, we present different tests our code was put through to benchmark its performance and limitations. 
We focus on tests applying to the modified gravity sector. For the results of standard N-body code tests for this code, 
see [42]. 



A. Spherical mass test 



In order to study how well the code reproduces the exact result for a spherical mass in DGP, we start out with 
a grid with Ng grid cells and Np particles. We arbitrarily assume a box size of Lbox = 200Mpc//i and set a = 1. 
All particles are moved into a spherical mass of radius R and uniform density, which then corresponds to a mass of 
M = PmoL^oK ~ ' lO'^^ Mq/H (assuming = 0.3). Given the density field assigned from the particle positions, we 
then use the relaxation solver to solve for ip, and measure the field values and radial acceleration throughout the box. 

For the purposes of this test, we can vary the two parameters /3 and rc in Eq. (2.6) independently. 1/(3 determines 
the strength of the additional force mediated by ip [e.g., Eq. (2.10)], which is attractive for positive f3 and repulsive for 
P < 0. For given /3, rc controls the Vainshtein radius, i.e. the scale where the self-interactions of the (p field become 
important. 

Fig. 1 shows the field solution (left panel) and relative deviation of the acceleration from the Newtonian value (right 
panel) as a function of distance from the center, r, for Ng = 512, i? = 10 grid cells, corresponding to 3.9Mpc//i, 
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Scale factor k [h/Mpc], scaled to L^_^^ = 256 Mpc/h 



FIG. 2: Left panel: Dimensionless RMS residual C — ^ (r^) and RMS of the field vs. a in cosmological simulations (top), 
and Vainshtein radius for a single particle in the simulations, in units of rceii (bottom). Right panel: Relative deviation of 
linearized DGP power spectra smoothed with Ts/rcdi = 1.0 from the unsmoothed simulations (solid lines, for box sizes from 
64Mpc//i [thick] to 400Mpc//i [thin]), and relative deviation of power spectra for lower resolution simulations [Ng — 128, left 
dashed curves, Ng = 256, right dashed curves). Also shown as vertical dotted lines are the maximum wavenumbers considered 
for each simulation, fcmax = fcNy/8. 



and rc/R = 800 « 3100Mpc//i. The agreement with the analytical solution (thick lines) is generally very good, 
except at large radii, where the periodic boundary conditions become important, and around r — R, where artifacts 
of interpolating a sphere onto a cubic grid become visible. At distances of order the size of the grid cells, the force 
resolution becomes worse leading to a growing scatter in the measured acceleration. Note that the analytical solution 
necessarily assumes an isolated mass, since the superposition principle cannot be used due to the non-linear ip field 
equation. For this reason, we added an arbitrary zero-point to the ip field in order to match the analytical solution in 
the left panel of Fig. 1. Of course, such a zero-point does not affect any observable quantity, such as the acceleration 
shown in the right panel. 

We tested the field solution for several different values of R, Ng, /3, and Vc- Given the caveats pointed out, we found 
very good convergence to the analytical solution. In particular, increasing R reduces the interpolation artifacts at 
r — R, while decreasing R reduces the impact of the periodic boundary conditions. Solutions for different values of (3 
are shown in both panels of Fig. 1. While the solid lines show the exact solution to the full equation, Eq. (2.13), the 
dotted lines show the linearized result, Eq. (2.10). Clearly, for the parameters chosen, the effect of the non-linearity 
is significant for r of a few R or less. 



B. Convergence and resolution tests 



In most test cases we studied, the density field is sufficiently smooth that the relaxation converges to the desired 
low tolerance. However, for density fields with considerable small-scale inhomogeneities, such as the cosmological 
density field at late times, we find that this level is not achievable. This is because the non-linearity of the ip equation 
is important as soon as the overdensity S becomes of order unity, which in the N-body simulation corresponds to ±1 
particle per grid cell. Equivalently, the Vainshtein radius for a single particle in our N-body simulation is of order the 
grid scale (see Fig. 2, bottom left). Thus, the ip field reacts much more strongly to noise due to the discreteness of 
particles than the Newtonian potential, leading to increased residual errors of the approximate solution. 

If we write Eq. (2.6) in code units (see Appendix A) as L{ip) = /, then for an approximate solution tpi, r = 
''ceii(-^('^0 ~ /) is the dimensionless residual, where rccii = ibox/-^g is the size of a grid cell. In the following, we 
use the RMS of the dimensionless residual, C = \/{r^ as a benchmark for the convergence of the field solution. 
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We performed tests with sine wave density fields of various wavelengths in order to determine what residual L is 
acceptable in our simulations (see Appendix B). In case of a pure sine wave, the non-linearity in Eq. (2.6) vanishes, 
and the exact non- linear solution is identical to the linearized solution. We found that residuals of £ < 10~^° are 
safe, as for residuals at that level, the errors in the solution are negligible compared to the unavoidable truncation 
errors from taking numerical derivatives on the grid. 

In order to reduce the noise in the density field in the cosmological simulations and improve the solution to the 
acceptable level of £, we increase the number of particles from (iVg/2)^ to A^^ , and smooth the density field entering 
the r.h.s. of Eq. (2.6) as described in Section III. With these steps, the solution converges with a residual £ < 10"^*^ 
for all box sizes (left panel of Fig. 2). Note that the dimensionless RMS residuals never exceed a fraction of ^ 10~^ 
of the RMS of the field solution. 

Since the force resolution in our simulations is in any case limited to scales above rceii, a smoothing on the grid scale 
is not expected to degrade the resolution of the simulations significantly. We checked this by running linearized DGP 
simulations given by Eq. (2.10) (Section II C) using the same smoothing of the r.h.s. of Eq. (2.6), and comparing the 
resulting power spectrum to that of linearized DGP simulations without smoothing, for the same initial conditions. 
The result is shown in the right panel of Fig. 2. As expected, the effect of smoothing increases towards larger k. The 
smoothing effect on the matter power spectrum is positive, since the smoothing removes power in the ip field on the 
smallest scales, and ip mediates a repulsive force. This leads to an enhancement in the matter power spectrum. 

The smoothing effect remains below 4% for k < fcmax = which we adopt as the maximum wave number 

considered for each box, where A;Ny = Ng/Lhox is the Nyquist frequency of the grid. From our studies with different 
smoothing radii = 0.8... 3.0, we found that the smoothing effects on full DGP simulations show a similar k- 
dependence as those of the linearized DGP simulations, but are smaller by a factor of ~ 0.6. This is understandable 
since the p field is suppressed in dense regions in the full simulations, reducing the effect of the smoothing. Taking 
into account this factor, we correct the power spectra measured in the full DGP simulations for the smoothing effects 
using the curves shown in Fig. 2 (right panel). Note that in any case these effects are at the level of few percent or 
less. 

In order to assess the effect of the finite grid resolution, we also performed simulations with Ng = 256 and Ng = 128, 
and Np — N^. Fig. 2 (right panel) also shows the deviation of the power spectrum in these low-resolution simulations 
from the Ng = 512 simulation with the same initial conditions. The vertical lines indicate fcmax = for each 

case, i.e. the maximum k we consider for each simulation box. Below /cNy/8, the deviations are less than 10% for 
either Ng = 256 or A^^ = 128. Note also that the resolution effects are independent of the type of simulation (DGP, 
linearized DGP, or GR). Hence, they cancel when measuring the deviation of P{k) between different simulation types, 
and we expect this deviation to be accurate to the percent level for the range in k we consider. 



C. Quasi-static approximation 

In solving for the brane bending mode Lp, we have assumed that all time-derivative terms in the full equation of 
motion of p (in terms of physical coordinates): 

Dp+^ [{Hpf - (V^V.^)2] = ^Sp, (4.1) 

can be neglected with respect to the spatial derivatives in Eq. (2.6). Up to second order, there are three time derivative 
terms which have been neglected: 

To = a'^dttP (4.2) 
07.2 

Ti = j^dttp-V^P (4.3) 

T2 = §(5*V,^)^ (4.4) 

Ti and T2 appear as Ti — T2 in the equations of motion. While we cannot rigorously prove that the quasi-static 
approximation holds, we can perform a consistency test by measuring the terms Eqs. (4.2)-(4.4) in the simulations, 
and checking whether they are indeed small. The time derivatives of p are calculated i n each grid cell using the two 
neighboring time steps. Wc then take the RMS of each term over the grid, || Ti \\= ^ (Tf ), and compare with the 
RMS of the spatial Laplacian, || </? || . Fig. 3 presents Tq and T\ — relative to the spatial derivatives for our largest 
and smallest box sizes, and shows that they are well under control. The relative magnitude of all time derivatives 
are of order a few 10~^ or less. This is in keeping with the expectation that dp/dt ^ H p, and k/a ^ H in our 
simulations. We conclude that the quasi-static approach is a self-consistent approximation. 
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FIG. 3: The time derivative terms Eqs. (4.2)-(4.4) of (y3 measured in the simulations relative to the spatial Laplacian, as a 
function of scale factor a for our largest and smallest simulation boxes. The time derivatives always remain 5 — 6 orders of 
magnitude below the spatial derivatives. 



V. COSMOLOGICAL SIMULATIONS 



We performed a suite of cosmological simulations for three different types of gravity: unmodified GR, corresponding 
to a smooth dark energy model with the same expansion history as DGP, referred to as QCDM; linearized DGP gravity, 
using the linearized ip equation of motion [Eq. (2.10)], corresponding to a time-dependent rescaling of Newton's 
constant; and full DGP, solving for the full non- linear (/? solution [Eq. (2.6)]. Comparing linearized DGP with full 
DGP allows us to study the Vainshtein effect in a cosmological setting. 

The cosmological parameters are those of the best-fit flat self-accelerating DGP model to WMAP 5yr data [48] and 
are summarized in Tab. I. We generated the initial conditions at a = 0.02 {z — 49) from a modified transfer function 
output of CAMB [49] for a flat ACDM model with f^A = 1 - ^m- The ACDM transfer function was corrected for 
small early-time modified gravity effects in the DGP model using the PPF approach, as detailed in Appendix C. 

The simulations were run with Ng = 512 grid cells on a side, and Np — 512^ particles, i.e. one particle per grid cell, 
to reduce the shot noise in the density field (see Section IV B). We performed six simulations each for four different 
box sizes, from ibox — 64Mpc//i up to Lbox — 400Mpc/ft, (Tab. II). On an 8-corc machine, the QCDM and linearized 



TABLE I: Cosmological parameters. 



TABLE II: Simulation type and number of runs per box size. 
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FIG. 4: Left panel: Power spectrum measured in the QCDM effective dark energy cosmology simulations (top), and halof it 
predictions. The lower panel shows the deviations from the halof it prediction separately for each simulation box size. Right 
panel: Halo mass function dn/d In M200 measured in the QCDM simulations (points) and the mass function prediction for 
QCDM after [51]. Also shown is the mass function for a ACDM model with the same primordial power spectrum. The lower 
panel shows the deviation from the prediction separately for each box size. 



DGP runs require ^lOh of computing time, while the full DGP simulations require about 300h. 

Fig. 4 (left panel) shows the combined power spectrum from all box sizes for the QCDM simulations, including 
bootstrap errors. Also shown is the non-linear power spectrum for QCDM, calculated from the linear power spectrum 
using the halof it procedure [50]. The power spectrum measured in each box is used up to fc,„ax = ^Ny/S (see 
Tab. II), and different boxes are combined weighting by volume. The lower panel of Fig. 4 (left) shows the power 
spectrum relative to the halof it prediction, measured separately for each box size. The power spectra measured in 
different boxes clearly agree within the errors, and the deviations from the halof it prediction are within the accuracy 
(^ 20%) expected from this fitting procedure (especially given the significant departures from ACDM of the simulated 
expansion history). 

Fig. 4 also shows the non- linear power spectrum for a flat ACDM cosmology fixed to the same rim, h, and primordial 
normalization, corresponding to a linear normalization of as = 0.66 today. Note that due to the different expansion 
history and the earlier onset of acceleration in DGP, the linear growth factor is suppressed by ~ 15% in QCDM 
relative to ACDM. The repulsive brane-bending mode in DGP suppresses growth further, leading to a suppression in 
the linear regime of ~ 21% relative to ACDM. In the following, we will always compare the DGP results with QCDM, 
so that the expansion history effects are taken out, and all deviations are strictly due to modifications of gravity. 

We have also measured the mass function of dark matter halos in our simulations. The spherical overdensity halo- 
finding and combination of different simulation boxes was done as described in [20]. As our grid resolution is the 
same as in the simulations analyzed in [20], while we have a factor of 8 more particles, we conservatively increase the 
minimum required number of particles from 800 to 6400. The corresponding mass thresholds are given in Tab. II. We 
define the halo mass M200 as the mass enclosed within a sphere of radius i?200 1 so that the average density within the 
sphere is 200 times the average matter density in the universe, pm — ^mPcrit- 

Fig. 4 (right panel) shows the mass function, nin_Af = dn/dlnM2oo, measured in QCDM simulations, in comparison 
with the prediction calculated from the linear QCDM power spectrum using the fitting formula from Tinker et al.[51]. 
The lower panel shows the relative deviations from the prediction, separately for each simulation box. Within the 
mass range accessible with our simulations, Af2oo ~ 10^^ — 10^^ M^/h, the agreement with the prediction and among 
different boxes is good. We also show the predicted mass function for a ACDM model with the same initial power 
spectrum in Fig. 4. Apparently, the number of halos above ~ 10^"* Mq/H is significantly reduced due to the suppressed 
growth in QCDM. 
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FIG. 5: Slices through a full DGP simulation with I/box = 64Mpc//i aX z = 0. Top left: density in logarithmic scale; top right: 
dynamical potential ^f; bottom left: brane-bending mode ip; bottom right: non-linear suppression of if field: y)NL = ip — ipL (see 
text). 

VI. RESULTS 

In order to get a visual impression of some of the physics in DGP N-body simulations, we show shoes through one 
simulation box at ^ = (Lbox = 64Mpc//i) in Fig. 5. The slices are 64 cells thick, and for each pixel we take the 
maximum absolute values of each quantity over the thickness of the slice, for better visibility. The density field (top 
left; in logarithmic scale) is difficult to distinguish visually from the QCDM result for the same run, as the effects 
on the power spectrum are at the ~ 15% level (Section VIA). The top right panel in Fig. 5 shows the dynamical 
potential exhibiting the potential wells of massive collapsed structures. The brane-bending mode ip is shown in 
the lower left panel. On linear scales, <^ ~ v'i is proportional to the potential (Section II C), but evidently does 
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FIG. 6: Left panel: Relative power spectrum deviation of full and linearized DGP simulations from GR (QCDM) at 2 = 0. 
The dashed points at high k extend the range from kNy/8 to fcNy/4 and are only intended to be indicative of the trend. The 
dotted line shows the deviation expected in linear theory, while the dashed line indicates the deviation predicted using linear 
theory as input for halof it. Right panel: Same as the left panel, but measured in the simulations at different redshifts. 

not follow the potential within deeper potential wells, making it appear smoother. To make this more clear, we show 
the difference of the if solution from the linear value [Eq. (2.10)], (^nl = if — fL, in the lower right panel of Fig. 5: f 
is suppressed (i^snl < 0) in overdense regions, showing the Vainshtein effect at work in a cosmological setting. Note 
that quite low-mass structures which are not conspicuous in the potential 4* already lead to a suppression of ip. See 
Section VIC for a discussion of the Vainshtein effect in dark matter halos. 



A. Power spectrum 

Fig. 6 (left panel) shows the relative deviation of the linearized and full DGP power spectra from the QCDM result 
at redshift 0. The error bars are obtained from the 6 separate runs using a bootstrap procedure. By comparing 
simulation runs with the same initial conditions, and then averaging the deviations, most of the cosmic variance 
cancels out and we are able to obtain significantly less scatter. Also shown is the predicted deviation in the linear 
power spectrum. Note first that in the self-accelerated DGP branch, the scalar field mediates a repulsive force, leading 
to a suppression of the growth of structure. Due to the scale-invariant modification of gravity in linearized quasi-static 
DGP [Eq. (2.10)], the predicted linear deviation is also scale-independent. On linear scales, k < O.l/i/Mpc, both 
linearized and full DGP simulations agree well with the linear prediction. 

Apparently, the full DGP result departs significantly from linearized DGP on quasi-linear to non-linear scales. The 
Vainshtein effect begins to operate wherever overdensities become of order unity, and suppresses the deviation from 
GR. Note that small effects of the non-linear ip equation can already be seen on quite large scales, k ^ 0.05ft./Mpc 
corresponding to r ~ 60Mpc//i, not far from the acoustic features in the matter power spectrum. While we do not 
expect dramatic effects on cosmological parameter constraints from BAO, a modeling of these non- linear effects will 
be necessary for precision constraints in the context of DGP and similar braneworld models. 

We also show the deviation of Pnl.dgp from Pnl,qcdm, where Pnl is calculated from the corresponding linear 
power spectrum using the standard halof it prescription, halof it describes the linearized DGP power spectrum 
reasonably well up to < 0.4/i/Mpc, in agreement with the findings of [52], while it follows neither the linearized 
nor full DGP at higher k. We have not explored whether phenomenological modifications of halof it allow for an 
improvement of the fit. 

The dashed points at high k extend the range in frequency up to fcNy/4 of our smallest box, 64Mpc//i. While 
resolution effects are still expected to cancel out to first order in this representation, these points only serve to 
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FIG. 7: Relative deviation of the halo mass function winA/ = dn/d\nM2oo measured in full and linearized DGP simulations 
from GR (QCDM) at z = (points). The dashed line shows the relative deviation of the predicted ACDM mass function (for 
the same Qm, h) from QCDM (see Fig. 4). 

indicate that the trends seen for full and linear DGP at lower k continue towards smaller scales. 

The right panel of Fig. 6 shows the evolution of the modified gravity effects on the power spectrum as function of 
redshift. On large scales k < 0.1/i/Mpc, the deviation is almost scale-free and evolves as predicted by linear theory. 
At earlier times, the density field is non-linear only on smaller scales. Hence, the Vainshtein effect becomes visible in 
the power spectrum only at higher k. However, it does affect the power spectrum deviation significantly already at 
z>l. 



B. Halo mass function 

The abundance of massive dark matter halos is a sensitive probe of the growth of structure in the Universe [SB- 
SB]. In particular, the number of the most massive halos which host galaxy clusters depends exponentially on the 
amplitude of the matter power spectrum. Fig. 7 shows the effect of the modification of gravity in DGP on the halo 
mass function, relative to the QCDM effective dark energy cosmology. We investigated the effect of the smoothing 
used in the full DGP simulations by comparing linearized DGP simulations with and without smoothing, as done for 
the power spectrum (Section IV B). Above our rather conservative mass threshold for each box, we only found a small 
effect of the smoothing on halo masses in the linearized DGP simulations. Since the ip field is in any case suppressed 
within halos due to the Vainshtein mechanism, we expect the smoothing effects to be even smaller for the full DGP 
simulations. 

As expected, the full DGP simulations are somewhat closer to QCDM than the linearized DGP case in Fig. 7. The 
suppression in the mass function is significant for M200 ^ 10^^ Mq//i, reaching more than 30% for massive cluster-size 
halos. The suppression is smaller for lower-mass halos. This is presumably because these halos formed earlier in 
cosmic history, when the modified gravity effects of DGP were significantly smaller. 

Note that when compared to ACDM, this suppression comes in addition to the larger suppression of the mass 
function due to the expansion history of DGP (see dashed line in Fig. 7). However, we expect the effect of the DGP 
modification of gravity to be fairly independent of the expansion history. In particular, we expect a similar effect on 
the halo mass function in generalized braneworld models which exhibit an expansion history close to ACDM [32, 38]. 
Hence, the mass function is expected to be a sensitive probe for braneworld modified gravity models, as is the case 
for ,f(R) gravity [20]. Note that for the normal branch of braneworld gravity, the sign of the effect in Fig. 7 will be 
reversed, leading to an enhancement of the number of massive halos. 
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FIG. 8: Left panel: Average halo profiles of 95, the dynamical potential ^, and the lensing potential $_ = ^fjv in the DGP 
simulations at z = 0, for halos between 10^'^ and 10^'' Mq/H. The lower panel shows the average profile of the matter overdensity 
6 for the same halos. Right panel: Average halo profiles of the radial gradient of ip, the DGP potential and Newtonian 
potential 'i'N, showing that gr — — Vr^' gN = — VrVl'iv in the inner regions of halos, i.e. GR is restored. Also shown is the 
brane bending mode calculated in the approach of [38], <^kw (see text), for the same halos as in the left panel. 

C. Brane bending mode and Vainshtein effect 

Fig. 8 (left panel) shows the average profiles of the brane-bending mode (p, the dynamical potential "if, and the 
lensing potential of dark matter halos in the full DGP simulations. The brane-bending mode has been scaled by 
3/3/2, so that it agrees with = "^n for the linearized solution. We only use our highest resolution simulation box 
for this measurement (Lbox = 64Mpc//i). The inner radius of the profiles is given by one grid cell. The profiles were 
measured by spherically averaging many lines of sight for each halo with 10^^ Mq/H < M200 < 10^'^ Mq/H, and then 
stacked, scaling each profile to the radius r2oo of the respective halo. This was done in order to reduce the significant 
scatter in the profiles. 

While the lensing potential $_ obeys the standard Poisson equation in DGP, the dynamical potential ^ receives a 
contribution from the brane bending mode [Eq. (2.8)]. Apparently, the if field flattens out towards the halo center, 
which is in qualitative agreement with the non-linear field solution for a spherical mass (see, e.g. Fig. 1). Note that 
(p turns away from the linear solution when the overdensity is of order a few, in agreement with the estimate in 
Section lie. 

The quantity which is actually observable however is the gradient of the field, shown in Fig. 8 (right panel), which 
gives essentially the deviation of the acceleration from the Newtonian acceleration. For each halo, we measure the 
radial gradient relative to the halo center. We again scaled V^V? by 3/3/2 to match the linearized solution to Wr'^N- 
From the halo exterior towards the interior, the gradient of ip grows initially, turns around at r ~ 0.8 r2oo, and shrinks 
for even smaller radii. As is clearly shown by comparing the dynamical potential ^E" with its Newtonian value, ^^at, any 
observable deviations from GR are indeed suppressed within massive halos. We find that for halos between 10^^ and 
10^^Mq//i, the deviation in the acceleration is around Ag/g^ ~ 0.2/(3/3) in the inner parts, while it is suppressed 
down to Ag/gN w 0.05/(3/3) for halos around 10^^Mq//i. This trend with mass is expected, given that cx M^/^ 
for the spherically symmetric solution. 

In order to model the behavior of the brane bending mode in the cosmological context, we can make use of the two 
limiting cases presented in Section II. One is the spherically symmetric solution (Section II D), which can be used in 
modeling a spherical collapse of dark matter halos in DGP [57]. Moreover, Khoury and Wyman [38] have based the 
approximate field solution in their N-body simulation on this limiting case. More realistically however, cosmic structure 
forms by collapsing into non-radially symmetric structures such as filaments. The plane wave density perturbation 
can serve as another, albeit rather extreme, test case. For such a perturbation, the non-linear self-coupling of (p 
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vanishes (Section TIE). Hence, we expect the Vainshtein effect to be typically weakened in a non-radially symmetric 
setting, compared to the spherically symmetric case. 

In their ansatz, Khoury and Wyman assumed that the spherically symmetric solution Eqs. (2.14)-(2.15) holds 
wherever the field becomes non-linear. Since in this solution, the gravitational constant is rescaled locally by 1 + A(r), 
the Poisson equation for ^ can be written in this ansatz as: 

V^* = AttG [1 + A{S)] a^pmd, (6.1) 

where A is now a function of the local overdensity. Substituting S ~ M/pmr^ in Eq. (2.16) [or Eq. (2.12)] we see 
again that (r*/r)'' ~ (-ffo?'c)^/3"^f^m<5, with an order unity coefficient. Khoury and Wyman [38] chose: 

(7)' = ^=^"™'^ ^ m = ^^{Vl+-e-l)/e. (6.2) 

Hence, in the DGP model we simulated, e « 0.3 d. Given the density field in our DGP simulations, we can solve 
Eq. (6.1) by Fourier transform as done in [38], subtract the Newtonian potential, and compare the brane bending 
mode (fKW from this ansatz with our numerical solution of the full ip equation. 

Fig. 8 (right panel) shows the result for stacked halo profiles of the radial gradient of vkw- While (fKW shows a 
qualitatively similar behavior to the full solution, the suppression due to the Vainshtein efi^ect appears to set in at 
significantly larger radii in (pKW than in the full solution. This is qualitatively in agreement with our expectation 
that non-radially symmetric configurations experience a weaker non-linear suppression. The fact that (^kw does not 
approach the linear solution at large radii is a result of the approximation Eq. (6.1), as was derived in the appendix of 
[38]. The precise tensorial structure of Eq. (2.6) restores (f to the linear solution at large distances r ry- It would 
be worth investigating whether this simplified (and computationally much less expensive) approach can be extended 
to recover the linear cp solution at large scales, an essential feature of the full brane-bending mode interactions. 

Although we did not compare our results with a full simulation based on the spherically symmetric approach of [38] , 
these results seem to indicate that this approximation overestimates the non-linear suppression of the brane-bending 
mode, which might afi:ect observables such as the power spectrum or halo mass function. We point out however that 
the crossover scale Tc in the braneworld-inspired model considered in [38] is an order of magnitude smaller than our 
Tc, and the non-linearity in Eqs. (6.1)-(6.2) only becomes efii'ective at much higher overdensities 6 = O(IOO), so that 
for the models studied in [38], the solution is possibly less sensitive to this approximation. 

VII. CONSTRAINTS ON THE SELF-ACCELERATING DGP MODEL 

Wc now briefly discuss the impact of our simulation results on cosmological constraints on the sclf-accclcrating 
DGP model without A. First, regarding, baryon acoustic oscillations, our results for the DGP power spectrum on 
quasi-linear scales show that non-linear effects are not significantly enhanced at the BAO scale. We estimate that any 
modified gravity corrections to the BAO scale are below the percent level, and hence well within the uncertainties of 
current BAO measurements [58, 59]. Including baryon acoustic oscillations will increase the power of the combined 
CMB and Supernova constraints on the self-accelerating DGP model [30] with non-zero curvature to 4 5ct [60] . 

Second, weak lensing measurements constrain the amplitude of the matter power spectrum today, which in this 
model is significantly smaller due to the suppression of growth by the earlier onset of acceleration and the repulsive force 
mediated by the brane-bending mode. The linear power spectrum normalization at z = in the models we simiilated 
is (78(QCDM) = 0.56 and C78(DGP) = 0.52, while for a ACDM model with the same primordial normalization, it 
is cTs (ACDM) = 0.66. We found that the non-linear matter power spectrum in the full DGP simulations is always 
below that of QCDM (up to fc < 6/?,/Mpc, Fig. 6). Hence, the suppression of the non- linear power spectrum in the 
self-accelerating DGP model corresponds to a reduction in the inferred linear normalization erg today of at least 0.1, 
with respect to a ACDM model with the same initial conditions. Such a deviation is disfavored by about 1.5 standard 
deviations by current weak lensing measurements [61]. 

The abundance of massive dark matter halos is probed by cluster surveys. We showed in Section VI B that the 
abundance of halos above 10^'^ Mq/H is suppressed by around 20% relative to QCDM, which itself only predicts half 
as many halos in that mass range as a ACDM model with the same primordial power. The latter again corresponds 
to a shift in erg by 0.1, which is constrained to more than 4cr by X-ray cluster measurements [56], when taking into 
account the systematic errors. 
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VIII. CONCLUSIONS 

The N-body simulations of DGP gravity we presented here show how the self-interactions of the brane bending 

mode Lp, which arc responsible for restoring General Relativity in dense environments, influence the formation of 
structure in the universe. In the self-accelerating DGP model we simulated, this scalar field mediates a repulsive 
force, effectively weakening gravity. We indeed find that the field solution turns away from the linearized solution 
whenever the matter overdensity 5 becomes of order a few, and that the repulsive force is suppressed by more than 
an order of magnitude compared to its linearized value in the center of massive halos. We also compared our full 
solution to that obtained in the approximate ansatz of [38]. While their ansatz agrees well with our results in the 
densest regions, it seems to overestimate the non-linear suppression in less dense regions, such as the outer regions and 
environments of halos. This is in line with our finding that the non-linear suppression of tp is weaker for non-spherically 
symmetric configurations. 

The non-linear matter power spectrum in the DGP simulations shows that the suppression due to the repulsive 
field is amplified on non-linear scales for the linearized DGP simulations. In the full DGP simulations, this suppression 
is smaller, and eventually turns around on Mpc scales. The deviations between linerized and full DGP power spectra 
are noticeable already on quasi-linear scales, k > 0.1 /i/Mpc. At the BAG scale, modified gravity effects on the power 
spectrum are at the percent-level. 

We found that the abundance of massive dark matter halos is significantly suppressed in the DGP simulations, 
compared to a standard gravity simulation with the same expansion history. Again, the non-linear suppression of 
the (f field alleviates this suppression, and hence has to be taken into account when using observations to constrain 
this type of modified gravity. The effect on the halo mass function is in fact large enough to make this an interesting 
observational probe of more general braneworld scenarios. 

Independently of the CMB constraints [30], our results on the non- linear structure formation strongly constrain 
the self-accelerating DGP model (without A). In the future, we plan to extend our simulations to the normal branch 
of DGP, and more general braneworld (-inspired) models [32, 34, 38]. The key part of the simulations, solving for 
the non-linear interactions of the brane-bending mode, is generic to all of these models, and hence the code will be 
readily generalized to these cases. We also plan to attempt a modeling of the modified gravity effects, in the context 
of the halo model for example. This model can then serve as a framework for cosmological constraints on braneworld 
gravity which take into account the non-linear mechanisms inherent to this model consistently. Simulations of both 
f{R) gravity and DGP have shown that most interesting phenomena appear on scales of 1 to tens of Mpc, and in the 
abundance and environments of clusters. Fortunately, these scales are well accessible to observations, e.g., through 
weak lensing, the Lyman-a forest, and cluster surveys, which can then be used as precision probes of gravity on 
cosmological scales. 
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APPENDIX A: CODE IMPLEMENTATION AND DISCRETIZATION 

This appendix describes details of the N-body code implementation, focusing on the relaxation solver for Eq. (2.6) 
as the main non-standard part of the code. The code is written in C++, uses OpenMP for parallelization of most 
time-critical operations, and employs FFTW [62] for the Fast Fourier Transforms. 
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1. Code units 

The comoving code units follow the convention used in [63, 64], where 

x = a 1 — , t=-, p = a— , * = ^, p = a-^ — , Al 

To to Vo Vo PmO 

ro = to = Hq^, Vo = ^, PmO = Pcrit.O, V'O = ^'o' (-'^2) 

In the above definitions, .Lbox is the comoving simulation box size in Mpc//i, Ng is the number of grid cells in 
each direction, Hq is the Hubble parameter today, Pcrit.o is the critical density today, and fim is the fraction of 
non-relativistic matter today relative to the critical density. 
Transformed to code units, the equation for (f becomes: 

Here, V acts with respect to code coordinates, fc = Tc/ro, c = I/vq, and S = dp/'p^ is the matter overdensity on 
the grid, determined from the particle positions. Note that we have not yet scaled ip to i/jq. After solving Eq. (A3), 
tp is added to the standard Newtonian potential, which is solved for using a Fast Fourier Transform, to obtain the 
dynamical potential in DGP: 

^ = ^N + ~- (A4) 

Z tpo 

2. Discretization of ip equation 

We employ standard second-order symmetric differences for the discretization of the second derivatives in Eq. (A3): 

VxVx<Pi,j,fe = {Vi+l,j,k + ^i-l,j,k - '^'fii,j,k) , (A5) 

^x^yVi,j,k = -^{'fii+l,j+l,k-Vi+l,j-l,k-^i-l,j+l,k + ^i-l,j-l,k), (A6) 

and correspondingly for derivatives with respect to y, z. Here, k stand for the grid indices, and the step size h = 1 
for the base grid and h = 2" for the grid of refinement level n. The finite differences are evaluated before squaring in 
calculating the non- linear terms in Eq. (A3). 

We have tried different discretizations, such as going to higher order in the finite differences, and solving for 
the deviation of (p from the solution of the linearized equation [Eq. (2.10)], instead of solving for Lp itself. The 
performance of those discretizations is comparable to or worse than the simple discretization Eqs. (A5) (A6). This 
is understandable, since going to higher order amounts to making the derivative operations less local in order to be 
sensitive to error modes of longer wavelength. However, in our multigrid relaxation scheme, this is already done 
efficiently by the coarser grids, so going to higher order in a single relaxation does not improve performance. Hence, 
we stay with the straightforward discretization, Eqs. (A5)-(A6). 



3. RelcLxation algorithm 

In general, a relaxation scheme operates by iteratively obtaining a better approximation to the solution (^('+^) given 
a previous guess, (p^^\ On the grid, the solution is updated successively for each cell {i,j,k): 

^IS^ - ^(^li. J' (A7) 

where F is the discretized field equation solved for tpij^h- In case of a linear field equation, F can be solved for ipij,k 
straightforwardly. However, in our case the field equation is non-linear in ip, and we instead solve for F iteratively 
via Newton's method. Writing the field equation Eq. (A3) as L(ip) = /, we need to solve: 



L{. ■ . , 'Pi,j,k, ■ ■ ■) — fi,j,k = 



(A8) 
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for ipij^k, where we have suppressed the dependences of L on neighboring grid cells. By expanding L{(pij^k) in a 
Taylor series around the current approximation, one step of Newton's method works by updating 'Pi'j^k with: 



^ <pw _ _ '-'J,- ' _ (A9) 



Our relaxation step is thus given by Eq. (A9), where dL/dtpi_j^k is determined from Eq. (A3) and the discretization 
Eqs. (A5) (A6). In practice, the order in which wc loop over grid cells in performing the relaxation Eq. (A9) is 
important, since each relaxation step depends on the neighboring grid cells. In particular, we need to take care when 
parallelizing the relaxation. We use a generalized red-black scheme (see, e.g., [42, 47]), by successively running over 
cells with [i+j + k) modulo 4 = n, n = 0, 1, 2, 3. This was done in order to break dependences due to neighboring cells 
within each set, in particular due to the mixed derivatives Eq. (A6), so that each set can be efficiently parallelized. 
We experimented with different ordering schemes and found that this choice performed best. 

Now, assume we have an approximate solution if to the field equation Eq. (A3) on the fine grid Vt^ , which differs 
from the true solution (p by the error e, = — e. Further, the residual is defined as r = L{if) — f. Since L{ip) = f 
by assumption, we obtain the following residual equation: 

L{ip) - L{ip -e)=r. (AlO) 

We expect that the relaxation on the fine grid has removed most small-scale error modes, so e mainly consists of 
longer wavelength modes. Thus, we will solve for e on the coarser grid fl^^, and then correct the approximate solution 
if for this error. On the coarse grid fP^, Eq. (AlO) reads: 

^2h^2h^^h^ _ L2h(^j2h^h _ g2h^ ^ j2h (i?^(<^ft) _ fh^ . (All) 

Here, the superscripts denote which grid a given quantity is defined on, and I^^ is the restriction operator mapping 
fields from the fine grid to the coarse grid. Note that Eq. (All) is of the same form as the cp equation on the coarse 
grid, L^''((p^'*) = p'^, except with a different right-hand side. Hence, Eq. (All) is solved for e^'* using the same 
algorithm as the original equation, but at one grid level higher. Once Eq. (All) is solved, we can correct ip^ for the 
error: 

jh^^2h^ (A12) 

where /2/i is the interpolation operator mapping the error from the coarse grid to the fine grid. 
In summary, the multigrid relaxation proceeds as follows: 

• Relax Eq. (A3) on the fine grid O'' Nj times. Compute the residual = L'^{(p'^) - 

• Solve Eq. (All) for the error e^'*: 

- Relax Eq. (All) on fl"^^ Ni times. Compute the residual (l.h.s.-r.h.s.) of Eq. (All). 

- Solve Eq. (All) for e'^^: 

- Correct e^'' using /f^e^'' [Eq. (A12)]. 

• Correct ip'' using ^^e^'^ [Eq. (A12)]. 

Thus, the relaxation proceeds as a nested loop going down to the coarsest grid, and then correcting the solutions on 
the successively finer grids. One such iteration through all grid levels is called a V- cycle. The coarsest grid we use 
has 4 cells on a side, which corresponds to a refinement of 2*" for a iVg = 512 base grid. 

For the interpolation operator Q^'* ^ fi'', we use standard bilinear interpolation (consistent with the cloud-in-cell 
scheme used for assigning densities and measuring accelerations on the grid). For the restriction — > fl'^^, we use 
full-weighting, the transpose of the bilinear interpolation operator. See the appendix in [42] for an explicit definition. 

In order to ensure convergence at all times, we ad opt a large value of Nj — 10, at the expense of computing time. 
Usually, one V-cycle reduces the residual C = i/ (r^) by 1-2 orders of magnitude. We stop the relaxation when either 
C < 10""'^^ is reached, or the convergence stalls at a certain level of £. Either of this usually happens within 4 V-cycles. 
The value of C for each time step is logged, allowing for a monitoring of the convergence status in the simulations. 
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FIG. 9: Truncation errors (red points) and residual errors in the ip field solution (green, blue points) for sine waves with 
k = fcNy/32 (left panel) and k = fcNy/16 (right panel). The parameters used for the field equation are shown (with rc in 
Mpc/h), as well as the RMS residuals €. of the final solution. The errors are shown as deviations of the acceleration measured 
at random points from the expectation, scaled to the RMS acceleration of the sine wave (see text). 



APPENDIX B: SINE WAVE TEST 

In this section, we use the fact that the full if solution is identical to the linearized solution for a plane wave 
perturbation (Section HE) to test our code for perturbations on different length scales. We consider a simple sine 
wave density perturbation, with fc-vector chosen to lie along the a;-axis: 

S(x) = — = Aainkx. (Bl) 
P 

In this case, the solution to the non-linear ip equation is identical to the linear solution: 

Hence, the exact acceleration is given by: 

- (1 + ^) M^) = (1 + ^) ^oskx. (B3) 

For a given level of residuals C, we now demand that the errors due to the approximate solution of the ip equation 
of motion are small compared to the unavoidable truncation errors which are encurred by taking the gradient of the 
potential on the grid. 

In the following, we set a = 1, = 0.3, and the number of grid cells is set to 256'^. The truncation errors are 
measured by assigning the analytical solution for the potential to each grid point, and comparing the exact acceleration 
to the one obtained by derivation and interpolation from the grid. The truncation errors for two sine waves of different 
wavelength are shown as deviations in the acceleration in Fig. 9 (red points). We scale the deviations to the RMS 
acceleration of the sine wave, (?rms = Sflm/i'^V^a k) A, in order to avoid far outliers obtained near the zeros of 
the sine wave when scaling to the exact solution at each point. For shorter wavelength modes, the truncation error 
increases and reaches order unity for waves on the Nyquist scale, as expected. 

In order to estimate the additional error in the acceleration made due to insufficient convergence of the solution, 
we compare the acceleration measured from the non-linear field solution coming out of the relaxation solver, to the 
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FIG. 10; Deviation of the DGP and QCDM matter transfer functions from ACDM at the initial redshift of the simulations, 
Zi = 49. The solid line shows the PPF calculation for DGP, while the blue dotted line is a simple fit used to correct the initial 
conditions (see text). 



corresponding solution of the linearized (p equation Eq. (2.10). We vary the parameters (3 and Vc which control the 
strength of the non-linearity in Eq. (2.6). In Fig. 9, we show the measured deviations in the acceleration from the 
linear solution for different values of P and rc, and in particular for the largest values that reached convergence. The 
final residuals of the (p solution in each case are also shown. 

Apparently, the errors in the acceleration are completely negligible for residuals C < lO"^**, consistent with results 
we found for other test cases such as the spherical mass. For sufficiently strong non-linearity and short wavelengths, 
residuals of order C ^ 10~^^ — 10~^^ are reached which lead to measurable deviations in the acceleration. However, 
for the residual values up to ~ 10"^*^ that we explored, the error due to the incomplete convergence is very small 
compared to the truncation error in all cases. Hence, we conclude that for residuals smaller than the conservative 
bound of 10~^°, the ip solution we obtain is sufficiently accurate in order not to bias the particle dynamics. 



APPENDIX C: CORRECTING INITIAL CONDITIONS FOR DGP 



Due to the additional term H/rc in the Friedmann equation, the self-accelerating DGP model starts to deviate 
from a ACDM expansion history at relatively high redshifts. In addition, there are modifications to the growth of 
horizon-scale modes. For this reason, even at the initial redshift of our simulations, Zi = 49, there are small departures 
in the matter transfer function in DGP as compared to a ACDM model with the same early-Universe parameters 
ilm = 1 - ^A, ^b, h, Us, Ag. 

In order to take these differences into account, we calculate the dark matter transfer function T{k) at Zi for DGP 
and QCDM, using the PPF approach described in [30, 65], and for the corresponding flat ACDM model. We do not 
include the effects of radiation on the transfer function, which is not necessary since we are only interested in the 
relative deviation of the DGP T{k) from ACDM. Fig. 10 shows the relative deviation in T{k, Zi) from ACDM for 
DGP and QCDM. Both QCDM and DGP transfer functions are slightly suppressed on all scales due to the effects of 
the earlier onset of acceleration. On super-horizon scales, the transfer function is further suppressed in DGP caused 
by the transition to 5D gravity on very large scales. In contrast, T{k) is less suppressed in QCDM because of the 
fluctuations in the effective dark energy. 

We correct the transfer function obtained from CAMB for these small differences at Zi by multiplying TcAMB(fc) 
with 1 -f /(fc), where f{k) is a simple arctan-t-const. function fit to the DGP curve shown in Fig. 10. Hence, these 
initial conditions are not quite correct for the QCDM simulations. However, the differences are very small on the 
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scales probed by our simulations {k > 0.01 h/Mpc). 
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